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We calculate the screening of the ion-ion potential due to electrons in the presence of a large 
background magnetic field, at densities of relevance to neutron star crusts. Using the standard 
approach to incorporate electron screening through the one-loop polarization function, we show 
that the magnetic field produces important corrections both at short and long distances. In extreme 
fields, realized in highly magnetized neutron stars called magnetars, electrons occupy only the lowest 
Landau levels in the relatively low density region of the crust. Here our results show that the 
screening length for Coulomb interactions between ions can be smaller than the inter-ion spacing. 
More interestingly, we find that the screening is anisotropic and the screened potential between two 
static charges exhibits long range Friedel oscillations parallel to the magnetic field. This long-range 
oscillatory behavior is likely to affect the lattice structure of ions, and can possibly create rod- like 
structures in the magnetar crusts. We also calculate the imaginary part of the electron polarization 
function which determines the spectrum of electron-hole excitations and plays a role in damping 
lattice phonon excitations. We demonstrate that even for modest magnetic fields this damping is 
highly anisotropic and will likely lead to anisotropic phonon heat transport in the outer neutron 
star crust. 
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I. INTRODUCTION 



Highly magnetized neutron stars called magnetars that feature extreme magnetic fields (B) — as large as 10 15 G 
on the surface and perhaps even larger fields inside — have been detected in recent years. Currently there are 
13 magnetar candidates which arc classified as anomalous x-ray pulsars (AXPs) or soft gamma repeaters (SGRs) 
(| http: / /www.physics.mcgill.ca/^pulsar/magnetar/main.html| ). Such large magnetic fields can strongly influence the 
state of matter in neutron stars, especially in the lower density region were the characteristic energy scales associated 
with the matter fields is small compared to the magnetic "perturbation" . Several decades ago, it was realized that 
even modest magnetic fields could influence the structure of atoms in the atmospheres of magnetized neutron stars 
and affect the spectral features of the emitted radiation [l| . At the surface, where electrons are localized around 
nuclei and the system consists of neutral atoms, large magnetic fields distort the atomic structure resulting in rod-like 
or cigar shapes elongated parallel to B. (For a recent review see 0.) It was realized quite early that these cigar 
shaped atoms would bind with each other along the magnetic field direction forming polymer chains, and that these 
chains would interact in the perpendicular plane, resulting in a condensed phase [3j. More recently the physics of 
these condensed phases for hydrogen has been analyzed and its implications for neutron atmospheres is discussed in 
Ref. H| . Here we will consider the effects of the magnetic fields at higher density, in the crusts of magnetars where 
the electrons are not localized in atomic states but rather form a degenerate Fermi system. This will be important 
in determining the structural and transport properties of the magnetar crust. Our finding will have implication for 
crust oscillations (excited in explosive phenomena such as giant flares) and for heat and electrical transport which are 
important ingredients in the study of magnetic and thermal evolution of magnetars 0, Q . 

In this work we will typically restrict the study to strong fields where only a few Landau levels are occupied and 
this will necessarily restrict us to relatively low density region in the outer crust of the neutron star, [p < 10 10 
g/cm 3 for B < 10 15 G.) The key microscopic quantity required to characterize the response in our approach is the 
polarization function, which is the Fourier transform of the density-density correlation function. The expression for 
the polarization function for a magnetized electron gas is calculated and is found in agreement with previous work Q ■ 
Our new findings in this work are related to identifying and understanding the implications of the electronic response 
in strong magnetic fields for neutron star structure and transport. The real part of the electron polarization function 
is related to the screening of the ion-ion potential in the crust. In an unmagnctized crust, screening is dominated 
by the usual Debye screening. In neutron stars where the typical magnetic field is small (B < 10 11 G), the electron 
screening length is larger than the inter-ion distance and this implies that the effects of electron screening are usually 
unimportant in determining the crystal structure and thermodynamics. In contrast, we find that for large fields when 
only one or a few Landau levels are filled, the screening length can be either much smaller or much larger than in the 
unmagnetized case depending on the density and magnetic field. At relatively low density the screening is enhanced 
resulting in a screening length that is generically smaller than the inter- ion distance. Our finding of enhanced screening 
for large B is not new, this was already noted in early work which explored some its implications for nuclear fusion 
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reactions [8[ . The qualitatively new finding of our work is the realization that the sharp Fermi surface in the electron 
momentum distribution along the B field will give rise to long-range oscillatory behavior called Friedel oscillations 
in the ion- ion potential — a well studied feature in the condensed matter context for B = 0. The period of this 
spatial oscillation Xp = 7r/(fe|), where k z j the Fermi momentum of the electrons in the z-direction, where we have 
assumed that B is parallel to z. This potential falls off slowly as 1/z on a scale of several inter-ion distances and 
consequently can be very important in determining the phase structure of the crust. We speculate that this will lead 
to the formation of an anisotropic crystal with rod-like structures akin to the polymer chains expected in the atomic 
regime. 

The imaginary part of the polarization is related to dissipation and gives rise to the usual Landau damping observed 
in Fermi systems. In the crust of a neutron star, excitations such as lattice phonons couple to the electrons, and 
their mean free paths are largely determined by the efficiency of the electronic Landau damping where the phonon 
decays by producing an electron-hole pair. We show that this decay rate has a non-trivial dependence on the angle 
between their propagation direction and the magnetic field because the electron-hole excitation for large B is highly 
anisotropic. This anisotropy leads to an anisotropy in heat transport properties of lattice phonons. We find that the 
lattice phonons can be scattered more easily if they are moving perpendicular to the field, compared to when they 
are moving parallel to it. 

The paper is organized as follows. In Section [TT] we derive the polarization function of electrons in the presence of 
a magnetic field. In the Sections IIHI and HVl we apply the result to two physical problems. In Section ITO1 we present 
numerical results for Friedel oscillations. In Section llVI we discuss the effect of anisotropic polarization on the decay 
rates of the lattice phonons as a function of the angle between the propagation direction and the magnetic field. 

We discuss the implications for the structure of the magnctar crust in Section [V] Details of calculations of the 
expression for the screened potential are given in Appendix [A] and Appendix [B] 



II. POLARIZATION OF ELECTRONS IN A MAGNETIC FIELD. 



The Lagrangian for the electron gas in an external electromagnetic field and at finite chemical potential /i is given 
by, 

L = d 3 xip(i]p- m e + n^ip , (1) 



where D„ = du — ieA^ and m e is the electron mass. We consider a time independent, uniform magnetic field B in 
the z direction and choose a gauge such that the only non-zero component of the gauge field A is A x = —A x = By. 
The spectrum of electron energy levels in an external magnetic field is well known and electrons occupy Landau levels 
with energies given by 



E m = Vfe) 2 + 2meB + m 2 e , (2) 

where m is the quantum number associated with the Landau level and is the momentum in the z direction. The 
number density of electrons is 

r> poo i>oo 

n * = 75-52 ( / dkZ h + 2 E / dkZ f^ ( 3 ) 
where f m = (exp( " E "^T A1 ) + is the Fermi distribution function. At T = 0, we obtain 

n e = g(v^^I + 2 £ W - < - 2meB W ^- m i-2meB) . (4) 

m>0 

We will often restrict to the low-density or large magnetic field limit, where only the lowest Landau level is occupied. 
This occurs for 2eB > /j 2 — m\. 

The density-density correlation function is defined as 

XB{x^yn = \{T{n{x)n{y)}), (5) 

where n(x) = tfv (x)ip(x) is the density operator. The Lagrangian is time independent for a static external field 
and even though the Lagrangian looks position dependent in a particular gauge, the system itself is translationally 
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invariant. (The change in the Hamiltonian due to a translation in x can be undone by a gauge transformation.) 
Therefore we can write the correlation function in momentum space as, 

XB (x*,y») = I ^i£e^-^n B ( w ,q) , (6) 

where q 11 = (w,q) = (u, q±, q z ). 

To evaluate lis , we insert a complete set of many-body cigenstates between the two density operators in Eq. [5] 
We ignore the anti-particle contribution to the polarization function because they arc well below the Fermi level and 
retain only particle-hole contributions. In this approximation we can write the polarization as, 

eB f°° dk z 

n B(9 M )=^E/ "af fmO--fn)W(E m ,E n ,m,n,k z m ,k*,q ± ) 



III . II ' 



r 1 1 ii 

1 uj - E n + E m + ie lu + E„ -E m - ie l |fc »- fc ™ + «" ' 



where, m refers to the (initial) energy eigenstatc with energy E m = y/ {k m ) 2 + 2meB + m\ and n to the (final) state 
with energy E„ = \J (k z n + q z ) 2 + 2neB + to 2 . The function W contains the matrix elements of the density operator 
between different electron eigenstates and is given by , 

W = W(E m ,E n ,m,n,k z m ,k^q x ) = [(E m E n + k z m k z n + mg)g ro -i, w -i(g x ).H^_ 1|W _ 1 (g±) 



+ (V2meBv / 2neB)H m _ hn -i(q±_)H* ln (q 1 _) ^ 
+ (V2meBV2neB)H m , n {q ± )H: n _ 1 , l _ 1 (q ± ) 
+ (E m E n + k z m k z n + m 2 e )H„ hn (q ± )H* mn (q ± )] , 

where, 

m>n [ e -*V4^(i=i^£2i) L^(&P/2)] , Q . 
Here, L h a (x) are the associated Laguerre polynomials, and I = l/\eB is the magnetic length. For m = or n = 0, 

Hm-l.n-l(X) = 0. 

The imaginary and real parts of the response function can be separated by following the pole prescription, 

= p (— \ 1 =)=Fi*S(u>-AE). (10) 



uj-AE±ie y us-AE 
This gives the imaginary part, 

9m[n B (g")]=- (-Tr)V / ^/ TO (1 - f n )W(E m , E n , m, n, k z m , k z n , g x ) 

2 ^ 2?r (11) 

[<5(w - E„ + E m ) - 6(w + E n - E m )] 1*.*=*,^+,* , 

We want to consider situations where the temperature is much less than the chemical potential \i. In this case, it is 
appropriate to approximate f m by #(/i — E m ) and 1 — /„ by 6{E n — fi). 
The real part of TIb is given by, 

eB « r°° dk z 
Ke[n B (<f)]=— / ^f m W(E m ,E n ,m,n,k z m ,k* n ,q±) 

p[ 1 1 ]l 

tO — E n + E m U) + E n — E m 

In particular, to calculate screening, we will be interested in the static response, u = 0, which is given by, 

eB f°° Hk z 
3te[n B (0,q)]=— £ / -^f m W(E m ,E n ,m,n,k z m ,k z n ,q ± ) 



m,n 



(13) 
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In the limit T = 0, 

3te[n B (0 )q )]= — 2 J] _/ ^W(£?m,^n ) m,n,A^ ) fc*,g x ) 



(14) 

1 



V(fc g ) 2 +^ + 2mei3)P[ — ] | fe g= fcW . 

We will now use these results lor calculating two important physical processes of relevance to magnetars: (i) the 
screening of electric charge by the electron gas and the resulting modification of the ion-ion potential in the crust 
( Section IIII[) . and (ii) the anisotropic spectrum of electron- hole excitations (Landau damping) and its consequences 
for the phonon mean free paths in the outer crust (Section [TV). 



III. SCREENING IN AN ELECTRON GAS 

The screened potential between two ions with charges Z\ and Z 2 separated by a displacement r in an electron gas 
is given by 

Knod(r) = — g(r) , 

f d 3 q e iqr 

where g(r) = 4vrr J — (15) 

gives the modification ol the bare interaction due to the polarizability ol the electron gas. The function F(q) = 
— e 2 9fie[IIoo(w = 0,q)] is related to the static polarization function of the electron gas, IIoo(0, q). In evaluating the 
integral in Eq.[15l typically it is assumed that lor distances r> 1/fc/, where kf is the Fermi momentum ol the electrons, 
we can replace IIoo(0, q) by IIoo(0, 0). In doing so one obtains the usual Debye screened potential in coordinate space 

^mod(r) = 1 . 2 6 exp(-r/Aj) , (16) 



or g{r) = exp(— r/Xo) where Ad = 1/ \/F(0) is the Debye screening length, and is simply related to the Debye mass 
tod, by Ad = 1/mc. F(0) can be calculated using the well known compressibility sum-rule @ and is given by 

F ^ = c2 k- (17) 

where n e is the equilibrium density and /i e is the chemical potential and e 2 = 47r/137. 

For a free electron gas F(0) = e 2 kfn/n 2 with /1 = (fc? + to 2 ) 1 / 2 . The ion density nj = 3/(47ra 3 ) — where a is the 

average inter-particle separation (lor a simple cubic crystal a is related to the lattice spacing L by (47r/3) 1 / 3 a = L) 
— and the electron density are related by electric charge neutrality which requires that n e = Zni, where Z is the 
charge ol each ion. Consequently, the electron Fermi momentum kf and the average inter- ion distance a are related 
by k f = (Zi^Znj) 1 ^ = 5.66 {Z 26 ) 1/3 a -1 where Z 26 = Z/26. Further, from Eq.[17]we have Ad = 10.37^/vJ/kf where 
Vf = kf/fi. The ratio ol the screening length to the average inter-ion distance is therefore given by 

Ad 1.83 — 

V(W73^' (18) 

implying that in the relativistic regime, where kf 3> m e and Vf ~ 1, Ad — 1-83 a. Since the interaction between 
the nearest ionic neighbors is essentially unscreened in the relativistic electron gas, screening is not considered to 
be important for determining the structure ol the neutron star crust. This situation is changed il the electrons are 
non-relativistic. For Z 2 6 = 1, Ad < a for Vf < 0.3. However at lower densities, where Vf <§; 1, orbital effects from the 
localization of electrons may be important and a simple screening approximation is not valid. 

Extending the Debye approximation to an electron gas in a magnetic field, setting ^(q) to its long- wavelength 
limit Fb(c{ = 0), we can deduce how the Debye screening is modified by the magnetic field. When only one Landau 
level is occupied, we find 

AD(2eB > /i 2 - to 2 ) = Jv*, (19) 

eve-B V 
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FIG. 1: (color online) The screening length Ad versus the inter-particle separation a for Z = 26 and eB = 10eB c . The dotted 
line corresponds to Ad = a. For a > 512(fm), the lowest Landau level is occupied, and for large regions Ad < a . As a 
decreases, more levels are occupied. 



which can be smaller or larger than a. Hence the screening length between ions relative to their separation can tuned 
arbitrarily by introducing a large magnetic field. To illustrate this we plot Xd versus a for eB = 10eB c in Fig. Q] 
In the low density regime where only one Landau level is filled, Ap < a for a > 1164(fm). With increasing density 
(smaller a), Xd becomes larger than a and continues to oscillate around the B = value as more Landau levels are 
occupied. The amplitude of these oscillations rapidly decreases with n where n is the number of Landau levels filled. 
As noted earlier, for B = the screening length is always larger than a when electrons are relativistic. 

At low temperatures the Fermi surface of the electrons is sharp and this gives rise to non-analyticities in F(q). 
Therefore it is not possible to approximate F(q) in Eq. [15] by F(0) even for r S> l/kf. The screened potential in 
position space depends on the polarization function for values of q other than q = 0, and has qualitative new features 
compared to the Debye screened result. This phenomenon is well known for the free electron gas (B = 0) 

Kapusta and Toimela [l(| employed the fully relativistic polarization function for electrons at one-loop order to 
calculate the most general form of the screened potential in a free electron gas. The central result of their calculation 
is that F(q) has branch cuts along q = |q| = ±2fc/ + ir\ for 77 > 0. Consequently at large distances, where r ^ Xd 
the screened potential exhibits characteristic oscillations called Friedel oscillations. It is well established @ that in 
the non-relativistic limit the screening function is given by 

4£ 2 cos (2fc/r) 

5Fricdci(NR) (r > X D ) = ( 4 + g) 2 (r/X D ) 2 ' ^ 

where £ = m D /2kj = (e 2 (i)/(Tr 2 k f ). At very low density, where £ > 1 Friedel oscillations are non-negligible and have 
been well studied in condensed matter physics context 0. The potential between two ions for r < A^ is dominated 
by Debye screening. In the region r~Ao a simple analytic form does not exist and the screening function needs to 
be calculated numerically. 

It is also shown in Ref. [Io| that in the extreme relativistic limit, when m e — > 0, Friedel oscillations are modified. 
At large distances the screening function is of the form 



r9 ^3/2 sin(2fc/r) 

.9Fricdcl(R) - (^t,) , ^ • 

{f / Xd) 



Unlike in the non-relativistic regime, here £ = 0.0047 1, and amplitude of these oscillations is greatly suppressed. 
Furthermore they decay as 1/r 3 instead of 1/r 2 in the non-relativistic limit. Consequently, Friedel oscillations can be 
ignored in the relativistic regime. 

In the presence of a magnetic field, the screening function Fb(q) is no longer rotationally symmetric and we expect 
the screened potential to be anisotropic. We now proceed to calculate V(r) in this case. The calculation of the 
polarization tensor when the densities are high and several Landau levels are occupied, is technically complicated. On 
the other hand, in this limit one expects the results to be very similar to the free electron gas. Therefore we restrict 
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ourselves to the large field limit, when only the lowest Landau level is filled. In cylindrical variables we write F B (q) 
as a function of q z = qcosd and q± = qs'mO, where 9 is the angle between the q and the magnetic field B. Note that 
Debye screening — which depends only on F B (0) — is completely isotropic. 
In cylindrical spatial coordinates the potential between two ions 

v ( \ ZlZ2e2 ( \ 
V mo d{p,z) = 4?rr 9\P, Z ), 



r 



oo 



a iq z 



where g{p,z) = - dq J _q J _J (q A _p) / dq z , (22) 

Wo J -co q 2 + F B {q±,q z ) 



with r = \J p 1 + z 2 , F B (q±,q z ) — — e 2 ils(0,q), and IlB(0,q) is the static polarizability of the electron gas in the 
presence of a magnetic field. 

Restricting to m = in the sum in Eq. 1141 we have, 



2^ ^J-oo 2tt ^ e ' E n -E (23) 



The sum over n in Eq. [53] runs over all non- negative integers, but the most important contribution for eB 3> p 2 and 
eB 3> q 2 is the n = term. (This is shown in Appendix [X]) Therefore for simplicity, we drop the n > terms, which 
allows us to calculate analytic expressions for F B (q), and the screened Coulomb interaction in certain limits. With 
these approximations, 

e ,n ,eB, 2 ,, „ m f k f „ , 1 



F B (q^q z )~(-) 2 C4h- ql/i2eB) I S dk 
n 1 J-k* 



^/(k z +q z ) 2 +m 2 -^(^) 2 +m 2 

(z4J 



((k z + q z )k z + m 2 e + y/(fc z ) 2 + m\y/{k* + q z ) 2 + m 2 ) 
2^(k z ) 2 + m 2 ^J(k z + q z ) 2 + m 2 



where kj = \/ p 2 — m 2 is the Fermi momentum of the electrons in the z direction for the m = level. If higher levels 
are occupied, each will have a different Fermi momentum. 

Now we consider two limiting cases. First, in the non-relativistic regime where m e ^> k z p we obtain for large £?, 

^,w:)-(f).--'"^i , t , +i ;_, t ,,, ) 

' (25) 
e N 2/ei? N _„2 ,2 /9 , , 2fc| + q z , m e 



= (lY(^±) e -^/2 l0 , 



-2k z f 



The Debye mass is given by 



lim/ B (a, 9 >(-) 2 SS, (26) 

q->0 7T Z K, 



and since n e = e_B / (27r 2 ) \J /i 2 — m 2 , the above expression (Eq. I26|) is consistent with Eq. [TTJ The key feature of the 
expression for F B (q±,q z ) in Eq. [25] is its non-analytic behavior as a function of q z . In the complex plane, Eq. 1251 has 
branch cuts along q z = ±2fc^ + irj, r\ > 0. Restricting q z to the real axis, one sees a kink at q z = ±2fc| (Fig. [2]). 
On the other hand the expression is analytic in q±. Based on this we can expect long range oscillations in g(z) with 
wavelength Tr/kj- in the z direction, but no long range features in the x — y plane. 
In the relativistic regime, where m e \i we find that 



tt v 2 ' J_ k z 2\k z \\k z + q z \ K -\k z \ + \k z + q z \> 

j (f) 2 ( eB )e - 9 i/2eB q± < fc , ^ 

The debye mass square is = lim q _,.o F B (q±, q z ) = (f) 2 (^). The number of electrons in these limits is 
n e = eB / {2ir 2 ) jj? — m 2 ~ (eB / (2ir 2 )) ^ and therefore m|) = e 2 dn e /dp e = e 2 (eB / (2ir 2 )) as expected. 
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q z (MeV) 



FIG. 2: (color online) Fb as a function of q z for qx = 0, eB — 100e_B c , and three values of p. As the electrons become more 



relativistic (ji increases) 
limit. 



the kink at 2k% becomes less sharp. The shoulder at k% becomes non analytic in the ultra-relativistic 



Prom Eq. [27] we note that in the ultra-relativistic limit Fb does not have any non-analytic behavior at q z = 2k On 
the other hand Fb is non-differentiable at q z = kj. This behavior is different from that of a free electron gas, where 
the non-analyticity occurs at ±2fc| in both the relativistic and non-relativistic limits. However, this is an artifact of 
the ultra-relativistic limit, and for any finite m e , Fb is differentiablc at q z = kj. We plot Fs{q z ) for three values of 
p in Fig. [5] We notice that for fixed p there are two distinct features in the plot versus q z . First, there is shoulder at 
kf which becomes sharper as m e / fx decreases. Second, there is a kink at 2fc| which becomes less prominent as m e /p 
decreases. In the ultra-relativistic limit the kink disappears, and the shoulder becomes a non-differentiable point. 
However, for any finite m e , the fact that the kink becomes less prominent suggests that the Friedel oscillations will 
become weaker as m e / p decreases. Hence, we will now focus on the non-relativistic limit. 

In the non-relativistic limit, we can obtain an analytic expression for g(p, z) valid for z 3> We find 

g(p, z ) = 9d(p,z) + 9f{p,z) , (28) 

where, 



g D (p,z) = e- m °^+^ (29) 
is the Debye screening formula which comes from the pole in the integrand in Eq. f2Uat q = imp, and 



. . v /7T^cos(2£^) ml(ir/4)p 
9f[P,z) = - '- 



l{2k z f y+ml{l/2)\n{Ak z z) (30) 

x K 1 (pJ(2k z f ) 2 + m 2 D (l/2) ]n(4k z f z)) 

arises from the branch cuts at q z = i2kj + irj. K is a modified Bessel function of the second kind. The derivation 
of Eq. [28] is given in Appendix [Bl For p = 0, we see that <7f(0, z) exhibits long-range Friedel oscillations in the z 
direction given by 

g F (0,z) cos(2fc^z) . (31) 

Eq. [28] is a good approximation for z > For /)40 and z —> 0, we expect that g(p, z) should — > 1 since the 

short range behavior of V(r) can not be modified by screening. For z — > the expression Eq. 1281 breaks down because 
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(2kj) 2 + m 2 D (l/2) ln(4kfz) < but in our derivation we used the fact that z > 1/kj and we can not trust Eg. [281 for 
z < anyway. For p = our derivation remains valid but Eq. [25] can not be used directly since Ki(ap) goes as 
l/(ap) as p — > 0. The limit of g as p — > 0, however, is finite because apKi(ap) tends to 1 as p —> 0. Therefore we 
define, 

COs(2fcfz) TO 2 n (7r/4) 

9[P ' ' 7T (2kf) 2 +m 2 D (l/2)ln(4kp) " 1 j 

For z = the value of 3(2;) is well approximated by the Debye screened value 

g(p,z = 0) = e- m ^ . (33) 



r(fm) 




FIG. 3: (color online) Plots of g as a function of p for z = 0, and as a function of z, for p — for eB = 100 e_B c , ^ = 0.7 
MeV. g shows an exponential decay as a function of p, while along the z axis it shows long range oscillations oc — cos(2ktz). 
Also shown is the Debye screening length Ad (vertical line, brown online) and ai = (3/(47rn/)) 1 '' 3 (vertical dot dashed line, 
orange online, assuming Z = 26). The solid lines correspond to the analytic approximations given in Eq. [32] (p = 0) and Eg. 1331 
(z = 0), while the points correspond to a numerical computation of the Fourier transform in Eg. 1151 The right panel shows the 
3-d plot of g(p,z) in the p, z plane with z and p. 

We have calculated g(p, z) for eB = 100 eB c and p = 0.7 MeV numerically. For these values, only the lowest 
Landau level is occupied and the Fermi momentum of the lowest Landau level is kj = 0.48 MeV. So the system 
is mildly rclativistic and the electron Fermi energy is large compared to atomic binding energy. The two relevant 
length scales for the screened potential are the Debye screening length Ad = 468 (fm), and the wavelength of Friedel 
oscillations Xp = ir/kj = eB / (2im e ) — 1296 (fm). In the left panel (Fig. [3]) we plot g(p,z) along the z and the p 
directions separately. This allows us to see the Friedel oscillations in the z direction clearly, and also compare the 
numerical results with the approximate analytic expressions. Along the z axis, Eq. [32] is a good approximation to the 
numerical result, and in the x — y plane, the exponential Debye formula describes the result well. A 3-dimensional 
plot of g(p, z) in the p, z plane is shown in the right panel and depicts that the Friedel oscillations are restricted to a 
narrow cylinder of radius r ~ 1/y/eB in the x-y direction. 

For large chosen field of 4.4 x 10 15 Gauss, which is near the upper end of what can be found in magnetar surfaces, the 
amplitude of Friedel oscillations can be very large as shown in Fig. [3] Nonetheless, it is important to note that even 
small amplitude oscillations can affect the lattice structure due to their long-range character. In our zero temperature 
treatment the Friedel oscillations are undamped in the z direction due to our assumption of a sharp Fermi surface. 
A finite temperature will smear the Fermi surface and this will result in exponential damping of Friedel oscillations 
by the factor exp (— z/£t), where £t = 2nvFe/T is the thermal correlation length [§] and VFe = kj/fi is the Fermi 
velocity. In neutron stars, the temperature is small and the thermal correlation length £t ^> a where a is the inter- 
ion distance. Consequently even for small amplitudes, Friedel oscillations can have important effects due to their 
long-range nature when the screening potential from different ions add coherently. 
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IV. HEAT TRANSPORT BY LONGITUDINAL LATTICE PHONONS 



Lattice phonons (lPhs) are space and time dependent vibrations of the ions, which can transport heat from one 
region of a solid to another. Typically both longitudinal and transverse modes contribute to heat conduction. In 
this analysis we restrict ourselves to only longitudinal modes since they primarily couple to electron-hole excitations. 
At low temperatures, when Umklapp process are suppressed [III, an d the occupation numbers of lattice modes is 
small enough that the non-linear terms in the lPh lagrangian are small, the dominant scattering of lPhs is with the 
electrons. We first re-derive the relation between the specific heat conductivity of lattice phonons and the imaginary 
part of the electronic polarization tensor, and recall results for scattering of lattice phonons in a free electron gas. 
Then we consider scattering in the presence of a magnetic field. 

At low temperatures the form of the interaction term between lattice phonons and electrons is well approximated 

by 

Cei = ~ [ (Pr^ipdiC , (34) 

J el J 

where l/f e i = (Ze 2 ni)/ (m 2 Dy /mini), with nj, the number density of ions, mi, the mass of ions, and itid, the debye 
screening mass [ll| . 

The rate for phonon decay into an electron-hole (e — h) pair can be calculated from the thermal width of the phonon 
due to absorption by electrons, and is given by the imaginary part of the self energy correction due to the electron 
loop. 

wr (^) = -4 9fmn(«") , (35) 

J el 

where T(q^) is the inverse lifetime l/r(q lJ ') of the lPh, and q^ — (w,q). For an on-shell photon w = ciq. From Eg. [331 
we obtain the mean free path Z/(g M ) = r(g /i )c/, 

kM = -( ), (36) 

where q is the lPh speed. 

The thermal conductivity contribution of the lattice phonons from kinetic theory is n = ^C v cik. C v for lattice 
phonons is simply 

C = . (37) 

The final expression for the specific heat conductivity for a typical phonon with uj = 3T is then, 

K = ) = ( ^ 2T2 fl 1 

15 v w3?7in(w,q) ; v 453mn(w,q) ; ' v ; 

For B = the polarization function (ignoring the antiparticle contribution) of an electron gas is well known [9l.[l2l|. 
Here, for convenience, we will approximate S$mII(<j^) by its value at T = 0. The corrections to this approximation 
are small if T/ fx is small, except at the kinematic boundaries. At T = the imaginary part of n(g M ) is, 

ZmW)]=-^6{qv f -\u\). (39) 

For phonon decay and conductivity in Eq. [38] we are interested in the imaginary part for uj/q = ci < Vf, which is 
given by 5to[II((7 m )] = — /j 2 q / (2n) . 

For B ^ 0, Smpflf^)] is given in Eq. [TTJ For a representative large value of the field in a magnetar we choose 
eB = 10eB c , and show results for the imaginary part of the polarization tensor, and the specific heat conductivity n. 
The lattice phonons have a typical energy u> ~ 3T, and we take T = 5keV for our calculation. The magnitude of the 
momentum is given by q = tu/c s , and the z and the _L components are q cos 9 and q sin 6 respectively. The speed of 
lattice phonons depends on the depth, but we take a representative value c s ~ 0.05. 

We show the results as a function of the mass density, in commonly used units p\2 corresponding to 10 12 gm/cc. To 
convert the mass density, which is dominated by ions, to the electron number density and hence the electron chemical 
potential, we take the atomic number of ions to be Z = 26 and their mass number A = 56 (Fe). In Fig. 2] we show 
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FIG. 4: (color online) A plot of Srm[Il(cj = 3T, q z = — cos 0, g_i_ = — sin 0)] as a function of the angle between the magnetic 
field and the lPh propagation direction, for three different densities pi2(10 12 gm/cc). For this plot eB = 10eB c , T = 5keV, 
and c s = 0.05. cos# = corresponds to a lPh traveling perpendicular to the magnetic field, and cos# = 1 corresponds to one 
traveling parallel to it. At the lowest density, only the lowest Landau level is occupied and scattering is kinematically allowed 
only at one specific angle. As more levels are occupied, more angles are kinematically allowed. 



Smpgf^)] as a function of cos 8 for a three different values of p 12 . In Fig.[5]we show the specific heat conductivity 
parallel (k z ) and perpendicular to the magnetic field, as a function of pi2- From Fig. Uwe see that at the lowest 
densities, the electrons occupy only the lowest Landau levels, and the response is highly anisotropic, and peaked 
around very specific values of 8, where energy conservation and momentum conservation along the z direction can be 
simultaneously satisfied. As the density of electrons increases, they occupy higher landau levels and overall response 
is obtained by summing contributions from various levels, and is non-zero for several values of 8. 

To calculate n± and k z in Fig.[5j we approximate the effective 3m[ns(w, q)] for scattering parallel to the magnetic 
field by the average between 8 = [0,7r/4], and the value perpendicular the magnetic field by the average between 
8 = [7r/4,7r/2]. We see that the lPh conductivity is suppressed perpendicular to the magnetic field but not parallel 
to it. This result is easy to see in the ultra-relativistic limit. Energy-momentum conservation in this limit requires 
\k z +qcos 6\ — \k z \ = c s q, which can occur only at cos 8 = c s <C ^j. Only at large densities where inter- level transitions 
become kinematically allowed, does the net response function become relatively isotropic. Formally, n z is very large 
in our calculation, but we have ignored the Umklapp processes which will be important in the z direction, because 
electron scattering is kinematically suppressed. One can include them by using k™ 6 * = k^k u / (n e z + k u ), where k% can 
be read from Fig. [51 and k u is the Umklapp contribution. In the perpendicular direction we expect electron scattering 
to be dominant, and it is alright to ignore the Umklapp contribution. 

We note that averaging the response function over angles is a rough method to include the effect of anisotropic 
scattering. A rigorous procedure will involve solving the transport equation using the anisotropic collision term, that 
can be deduced from the expression for QmTls^)]. We leave such a calculation for future work. 



V. CONCLUSIONS 



By calculating the one-loop electron-hole polarization function in a strong magnetic field we have found that the 
screening of the ion-ion potential is significantly different from those obtained in earlier studies where the Debye 
approximation was assumed. We show for the first time that the Fermi surface of electrons in the direction parallel 
to the magnetic field leads to Friedel oscillations in the ion-ion potential. These oscillations are large and could 
fundamentally change the structure of the ionic solid at large field and relatively low density when only the lowest 
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FIG. 5: (color online) Plot of k as a function of density. k z is the specific heat conductivity parallel to the magnetic field, 
shown by the dashed (red) line. This should be considered as an upper bound for k z because our calculation ignores Umklapp 
processes. k± is the conductivity perpendicular the field shown by the solid line (green). The dotted line (black) is the 
conductivity for B = at the same density. 



Landau levels are occupied. For typical magnetar field strengths of order 10 15 G, matter up to densities of the order 
10 10 gm/cc will be affected by our finding here. To evaluate how Friedel oscillations and modifications to the average 
screening length affects the structure and melting properties of the solid we will need to perform either Path Integral 
Monte Carlo or classical Monte Carlo simulations because the ion-ion interaction is highly non-perturbative. This is 
beyond the scope of this work and will be investigated separately. Here, we present plausible implications based on 
qualitative arguments to explore how our findings will modify the ionic structure. 

Considering the parameters used in Fig. [3l eB — 100eB c and fi = 0.7(MeV) and assuming Z = 26, the average 
number density of ions in the system is nj = 3.2 x 10 -9 (fm) -3 . If one assumes that the ions form a regular bcc 
lattice, then the separation between nearest ions is (y/3/2) x (2/n,/) 1 / 3 = 743(fm). This is comparable to the screening 
length Ad, and smaller than the wavelength for Friedel oscillations A^ . If significant screening in the x — y plane 
prevents the formation of a regular bcc structure, it may be favorable to form an anisotropic crystal structure. One 
possibility is that the ions fit into the troughs in the potential formed by superposition of the V(z) separated in the 
z axis hy Xp. To maintain neutrality, the ions should arrange themselves in the x ~ y plane (in a regular or irregular 
structure) with an average separation aj_ = (nniXF)^ 1 / 2 = 278(fm). 

If indeed the ions in the z direction arrange themselves with separation Xp, then the following picture arises for the 
structure of the ion-electron system with changing density, or as we move deeper into the neutron star crust. A single 
chain of ions arranged along the z axis will look like a rod. Xf oc eB /nj, and for fixed eB decreases linearly with 
increasing density. The transverse separation between the rods, a± ex (eB) 1 / 2 , on the other hand remains a constant. 
At the surface, it is known that the elongated atoms form chains parallel to the B field. At low densities (near the 
surface of the crust), ions will form a plasma of rods with a small charge per unit length (since Xf is large). These 
rods will interact with each other by a two dimensional Dcbye screened Coulomb interaction. As we move deeper into 
the star, the charge per unit length will keep on increasing, and once Xf becomes comparable to aj_, the system will 
look more isotropic. This picture will break down when higher Landau levels start being occupied. It is interesting 
to note that for Z = 26, a± ~ \ F implies that /x 2 > 2eB, meaning that the next Landau level is occupied. Hence, 
the point where the screened potential becomes more isotropic coincides with the point where the lattice structure 
becomes more isotropic. Hence, in our picture, the region where the electrons occupy only the lowest Landau level 
interpolates smoothly between the surface and the inner crust of the neutron star. 

We have also showed that the heat conduction due lattice phonons can become anisotropic because their damping 
rate due to electron-hole excitations is anisotropic. It is well known that the electronic heat conductivity perpendicular 
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to the magnetic field is suppressed [13[ , and this in turn results in a temperature anisotropy at the surface of the neutron 
star [f|. Our finding that the heat conductivity due to phonons in the direction perpendicular to the magnetic field 
is also suppressed, is new. This effect may modify the temperature anisotropy in the outer regions of the magnetars 
during their early thermal evolution. 
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Appendix A: Dominance of the Lowest Landau level for eB 3> fJ?,q 2 
To see that the n — term dominates in the sum in Eq. 1231 we look at the integrand, 



W(E o ,E n ,0,n,k^,k^,q z ) _ (E Q E n + k z n k z n + m 2 )e-i 2l2 / 2 (q 2 l 2 /2) n /n! 1 



E m — E„ 2E E n E n — E„ 



(Al) 



where we have used Lq(x) = 1. Now, E m = J (k m ) 2 + m 2 ~ /j. For n = 0, E n = J (k z + q±) 2 + m 2 ~ p while for 
n > 0, E n ~ \j2neB. Therefore, for n = 0, 

W{E O7 E n ,0 7 n,k m ,k n ,q z ) ^ e - q ^/ 2 J_ ^ 
E n — E m 2[i 

while for n > 0, 

W(Eq, E n ,0, n, k m , k n ,q z ) e-i 2l2 / 2 (q 2 l 2 /2) n /n\ e -^ l2 / 2 (q 2 l 2 /2) n /n\ 



E n - E m 2(E n - E m ) 2\/2e~B 



(A3) 



Thus, we can conclude that the processes involving excitation to the level n is suppressed compared to the m = 0, 



n 



transition by a factor \/2enB)){q 2 / {2eB)) n (l/n\) which is small if eB > fi 2 , eB > q 

Appendix B: Derivation of g(r) in the non-relativistic limit 

To proceed with the derivation in the non-relativistic limit, we start form the expression, 



2 



= vV + 

where 



vV+Zf, j, v r,z exp(ig^) 
9iP,z = / dq±q±Jo{q±p) / dq 

k Jo J-oo [q z r + q± + F B {q±,q z ) 

/ dq_ L q ± J {q_ L p)I{q_ L ,z) , 
Jo 



(Bl) 



J-oo {q z r + (q±r + F B {q±,q~) 

We calculate the integral I(q±,z) by using contour integration. We assume z > and close the contour in the 
upper half of the complex plane. We have to deform the contour to go around the branch cuts ±2fcy + irj, rj > 
0. We call these two parts of the closed contour, C\ (q± = 2kf + irj) and C2 (q± = —2kf + irj). This gives 

We treat these contributions one by one. 

First, we look at the pole contribution which we will call the Debye part go- The pole is where 

= qi + (q z ) 2 + F B (q ± ,q z ) 

■e-,2,eB, 2 p ,„ ,2kf + q z ,m e (B3) 



tt> K 2 ' 2fc / - 2 q z 
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If we assume that the dominant contributions comes from q z <C 2k f, then we want 



(B4) 



= ql + ( q z ) 2 +m 2 D e~^' 2 



This gives, q z = iyjqj_+ m 2 D exp (—ql/(2eB)), where we consider the pole in the upper half complex plane. Fur- 
thermore if the values of q± that contribute strongly to the dqj_ integral are q± <C v / 2e£?, we can approximate 
exp {-q\/{2eB)) by 1 and 



exp (iq z z) 
lD{q_L,z) = 2ttzRcs( - ^ - ^ 2 ) 



(q z 



H+ m D 



2tt e -Wi±+ m l 



(B5) 



D 



This is exactly what we obtain if we replace F(q±, q z ) by m 2 D in the first place. Hence, this integral is the same as 
what we obtain if we assume that there is a simple Debye like screening, and therefore gives, 



/ID 



(p,z) = e~ mD V^+ z 



This is the answer if q± <C kf, \feB dominates the q± integral, and if (~ ) 2 (^j-)(m e /kf) <§C k 2 . 



Now we consider the integrals over the contours C\ and C2, 



dq 



o(iq z z) 



dq z 



Cl (q z ) 2 +ql+F B (q ± ,q*) J C2 (q*) 2 + q 2 ± + F B (q ± ,q*) 



da z 

Cl (q z ) 2 + ql + F B (q±,q z ) 



dr\ 2ie {2tk i-^ i 



(B6) 



ql + (2k f + t V ) 2 + (e/n) 2 (eB/2)(m/(2k f + ir,)) exp(-^/(2e J B))(l/2)(ln((4fc / + i V ) 2 /r/ 2 ) - in) 

1 

" q 2 ± + {2k f + i V ) 2 + (e/n) 2 (eB/2)(m/(2kf + i V )) exp(- g 2 _/(2eB))(l/2)(ln((4fc / + ir]) 2 /-!] 2 ) + in) 



ffte 



dr) 2ie {2ik '~ 7l)z ( K 



-1 



(B7) 



q\ + {2k f + i V ) 2 + m 2 D (k f /(2k f + in)) exp(- 9 2 /(2eB))(l/2)(ln((4fc / + ir,) 2 /ri 2 ) - in) 

1 

+ ql + (2k f + i V ) 2 + m 2 D (kf/{2k f + ir])) exp(-g 2 _/(2eB))(l/2)(ln((4fc / + ntf/ri 2 ) + in) 

Because of the cxp(— -qz) in the integration over 77, for z 3> 1/&/ the integral is dominated by r\ ~ 1/z <C kf. Therefore 
we can ignore r\ whenever it is added to kf. This simplifies the integral somewhat, 



lF(q±,z) ~ 5ie 



dr) 



2e (2ik f -, 1 )z m 2 D exp(-^/(2eB))(7r/4) 



(ql + (2k f ) 2 + m 2 D eM-ql /(2eB))(l/2) ln(4fc / /r ? )) 2 + ( m 2 D exp(-<£/(2eB))(7r/4))s 



(B8) 



To obtain an analytic form, we need one more simplification. We replace 77 in \n(4kf/r)) by the value where we expect 
the integral to dominate, namely rj ~ 1/z. The integral over 77 can then be evaluated simply, and 



lF(q±,z) 



2cos(2fc / z)m|,exp(-g 2 _/(2e J B))(7r/4) 



rl 

~z (ql + (2k f ) 2 + m 2 D exp(-g 2 /(2eB))(l/2) ln^fc/z)) 2 + (to 2 , exp^/peS))^))* 



(B9) 
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The Friedcl contribution to g is then, 

I r°° 

9f(p,z) = -vV + z 2 - I dqj_q±J{qj_p)I F (qj_,z) 
* Jo 

-1 2cos(2fc f z) f°° , 

= -^p 2 +z 2 — — / dq ± q ± J(q±p) BIO 



m 2 3 exp(- g 2 _/( 2e i3))^/4) 



+ {2k f Y + m 2 , exp(-£ /(2eB))(l/2) ln(4fc /2 r))a + (m 2 D exp(-gi/(2eB))(7r/4))a 

Once again assuming q± -C VeS we ignore the exp(— g^/(2e_B)) term as before. Finally, assuming kf 3> mjj we can 
drop the (mjjir/A) term in the denominator, and 

, \ r-o 5rl2coB(2fc/z) / ,00 J _. . mUn/A) 

gF{p , z) ~-^r^-——J o ^^^p) ((?i + (2fc/)2 + myi/2)ln(4 , /z))2 

(Bll) 



\J p 2 + z 2 cos(2kfz) m 2 D {-K/4)p 

* ~ z V / (2k f ) 2 + m 2 D (l/2)ln(ik f z) 



K, (p v /(2A; / ) 2 +m 2 ,(l/2)ln(4fc / z)) 



These manipulations work only if z 3> 1/fc/, (long distances in z) and p <C 1/fc/- The assumption p <C is not 
a very essential one because for p > 1/fc/, gp(p,z) decreases very rapidly and this is captured by the expression in 
Eq. IB111 because of the modified Bessel function, K\ . 

Thus, the final answer for large z is, 



g(p,z) = g D (p, z) + g F (p, z) 

— e -m D -\/p 2 +z 2 



J p 2 + z 2 cos(2fcfz) mn(7r/4)p . /. , .„ , , , " 

^ V /^nr ± I,,, =X 1 (p V /(2fc / ) 2 +m 2 D (l/2)ln(4fc / z)). 

2 V(2fc/) 2 +Wc( 1 /2)ln(4fc / z) v 



(B12) 
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